In [ ]:

    

#! /usr/bin/env python
# -*- coding: UTF-8 -*-
import numpy as np

import scipy.io as io
import scipy.stats as st
import scipy.linalg as la

import matplotlib.pyplot as plt
%matplotlib inline

import warnings
warnings.filterwarnings("ignore")


    

In [ ]:

    

mat = io.loadmat( "./Dog_1_interictal_segment_0001.mat" )
ndat = mat[ "interictal_segment_1" ]

l = ndat['data_length_sec'][0,0][0,0]
f = ndat['sampling_frequency'][0,0][0,0]
d = np.array( ndat.item(0)[0] )


    

In [ ]:

    

def prin1( d, l, f, M = 5, B = 2.5, N = 5 ) :
    L = N // B
## The time by which the directions are computed
    time = np.arange( 0, l + 1, B )
## The index in the original data (corrected for sampling frequency)
    index = np.floor( f * time )
## A 3D array to hold the principal directions
    pcv = np.empty( ( len( index ) - L, M, d.shape[0] ), dtype = np.float )
## A matrix of variances of the top M principal directions
    lam = np.empty( ( len( index ) - L, M ), dtype = np.float )
## A matrix of variances of the data
    var = np.empty( ( len( index ) - L, ), dtype = np.float )
    for s,t,i in zip( index[:-L], index[L:], range( lam.shape[ 0  ] ) ) :
## Compute the ususal correlation coefficient: its poor stabilty under outliers
##  plays to our advantage, since we do want amplified correlation in cases of
##  extremely volatile events (or onslaught of such).
        rho = np.corrcoef( d[:,s:t] )
        # rho, _ = st.spearmanr( d[:,s:t], axis = 1 )
## Since the matrix is symmteric, ordinary EVD suffices for the PCA
        e, v = la.eig( rho, right = True, check_finite = False )
## Since Numpy does not guarantee the descending order of the eigenvalues
##  sort them manually
        desc = np.argsort( e )[-M:][::-1]
        lam[i] = e[desc]
        var[i] = np.sum( e )
## Keep the eigenvectors corresponding to top eigenvalues 
        pcv[i] = np.transpose( v[ :, desc ] * np.sign( v[ 0, desc ] ) )
## Compute the geometric persistence of each principal direction:
##  the persistence is measured by the cosine between the previous
##  and the current principal direction
    align = [ np.diag( np.dot( pcv[:,k], pcv[:,k].T ), 1 )
        for k in range( M ) ]
## Return the whole bundle: time, eigen-vales and -vectors. Also give the 
##  cosines and the total variance.
    return time[L:], lam, pcv, align, var


    

In [ ]:

    

time, lam, pcv, per, var = prin1( d, l, f, M = 8, N = 20, B = 5.0 )

align1 = np.diag( np.dot( pcv[:,0], pcv[:,0].T ), 1 )
align2 = np.diag( np.dot( pcv[:,-1], pcv[:,-1].T ), 1 )


    

In [ ]:

    

import matplotlib.cm as cm
## Plot the variances of the top directions
plt.figure( figsize = (16,9))
colors = cm.rainbow( np.linspace(0, 1, lam.shape[ 1 ] ) )
for k in range( lam.shape[ 1 ] ):
    plt.plot( time, lam[:,k] / var, '-', color = colors[k], linewidth = 1 )

## Compare the largest and the tinyest directions and their geometric persistence
plt.figure( figsize = (16,9))
plt.subplot(211)
plt.plot( time, lam[:,0], '-k' )
plt.plot( time, lam[:,-1], '-r' )
plt.subplot(212)
plt.plot( time[1:], per[ 0 ], '-^k' )
plt.plot( time[1:], per[ -1 ], '-^r' )

## Show the directinal uniformity factor (the more uniformly the variance is
##  spread among the principle directions, the less structure in the multivariate
##  data there is).
plt.figure( figsize = ( 16, 9 ) )
plt.plot( time, lam[:,0] / lam[:,-1], '-k' )

plt.figure( figsize = ( 16, 9 ) )
ax = plt.subplot( 111 )
ax.set_color_cycle( cm.rainbow( np.linspace(0, 1, d.shape[ 0 ] ) ) )
original_time = np.arange( 0, d.shape[1] ) / f
for i in range( d.shape[0] ):
    ax.plot( original_time, d[i,:], '-' )


    

In [ ]:

    

import numpy.fft as fft


    

In [ ]:

    

freq = fft.fft( d[0:1], axis = 1 )


    

In [ ]:

    

sig = fft.ifft( np.log( np.abs( freq[ 0,: ] ) ) )


    

In [ ]:

    

plt.plot( np.reasig[:-1000] )


    

In [ ]:

    

plt.figure( figsize = ( 16, 9 ) )
plt.plot( np.log( np.abs( freq[ 0,: ] ) ) )


    

In [ ]: